#R code for robust regression models in the Supplemental Appendix, Transferring Power to Maintain Control#
rm(list = ls())
install.packages("lme4")
install.packages("influence.ME")
install.packages("haven")
install.packages("robustlmm")
install.packages("coeftest")
library(lme4)
library(lmerTest)
library(interplot)
library(haven) # for stata files
library("dplyr")
library(influence.ME)
library(stargazer)
library(robustlmm)
library(ggplot2)
library(directlabels)
library(RColorBrewer)
library(colorRamps)
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))
library(ggplot2)
library(directlabels)
library(RColorBrewer)
library(colorRamps)
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))
library(lme4)
library(lmerTest)
library(interplot)
library("dplyr")
library(influence.ME)
library(robustlmm)
library(stargazer)
library(quantreg)
options(scipen=999)

setwd("C:/Users/")  #you can set this to your own wd.

#load("Replication_data.RData")
l_data <- read_dta("Decentralization_WorldPolitics_Final_June_2023.dta")


#Model 1#
out1<-lmer(raiapproved ~ langdif_lg + area_lg + logpopmil_lg + jollyeuind_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out1)

out1rr<-rlmer(raiapproved ~ langdif_lg + area_lg + logpopmil_lg + jollyeuind_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out1rr)


##### To get p-values
##### See: https://stats.stackexchange.com/
##### questions/430490/obtaining-p-values-in-a-robustlmm-mixed-model-via-satterthwaite-approximated-dfs

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(out1)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(out1rr))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values


#Model 2#
out2<-lmer(raiapproved ~ econdif_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out2)

out2rr<-rlmer(raiapproved ~ econdif_lg +
                raiapproved_lg  +
                (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out2rr)


##### To get p-values
##### See: https://stats.stackexchange.com/
##### questions/430490/obtaining-p-values-in-a-robustlmm-mixed-model-via-satterthwaite-approximated-dfs

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(out2)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(out2rr))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values



l_m3data <- l_data[which(l_data$subtiminginbet_lg==1), ] 


#M3
out3<-lmer(raiapproved ~ coalgovtvote_lg + subnatcoalvote_lg + natcoalXsubnatcoal_lg +
                raiapproved_lg + 
                (1 | regionno ) + 	(1 | countryno), data=l_m3data,  REML=FALSE)

summary(out3)

out3_rr<-rlmer(raiapproved ~ coalgovtvote_lg + subnatcoalvote_lg + natcoalXsubnatcoal_lg +
                    raiapproved_lg + 
                    (1 | regionno ) + 	(1 | countryno), data=l_m3data,  REML=FALSE)
summary(out3_rr)


##### To get p-values
##### See: https://stats.stackexchange.com/
##### questions/430490/obtaining-p-values-in-a-robustlmm-mixed-model-via-satterthwaite-approximated-dfs

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(out3)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(out3_rr))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values

#Model4

out4<-lmer(raiapproved ~ ethnicthreat_lg + gvuln_lg + gvuln_lgXethreat_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out4)

out4_rr<-rlmer(raiapproved ~ ethnicthreat_lg + gvuln_lg + gvuln_lgXethreat_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out4_rr)


##### To get p-values
##### See: https://stats.stackexchange.com/
##### questions/430490/obtaining-p-values-in-a-robustlmm-mixed-model-via-satterthwaite-approximated-dfs

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(out4)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(out4_rr))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values


#M5

out5woEU<-lmer(raiapproved ~ ethnicthreat_lg + gvuln_lg + gvuln_lgXethreat_lg + langdif_lg + area_lg + logpopmil_lg +
             raiapproved_lg  +
             (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out5woEU)

out5woEU_rr<-rlmer(raiapproved ~ ethnicthreat_lg + gvuln_lg + gvuln_lgXethreat_lg + langdif_lg + area_lg + logpopmil_lg +
                     raiapproved_lg  +
                     (1 | regionno ) + 	(1 | countryno), data=l_data,  REML=FALSE)

summary(out5woEU_rr)


##### To get p-values
##### See: https://stats.stackexchange.com/
##### questions/430490/obtaining-p-values-in-a-robustlmm-mixed-model-via-satterthwaite-approximated-dfs

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(out5woEU)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(out5woEU_rr))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values

